###############################################################################################
## Replication R Code for
## No Reservations: International Order and Demand for the Renminbi as a Reserve Currency 
## Authors: Steven Liao, Daniel E. McDowell
###############################################################################################

###############################################################################################
## Code date: August 4, 2015
## R version 3.1.3 (2015-03-09)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Please contact Steven Liao <stevenliao@princeton.edu> for questions
################################################################################################

############
## Setup
############
# clear all objects in workspace
rm(list = ls())

# load packages
# sltool is a package in development by Steven Liao, 
# it includes various functions used in this study
# to install please follow instructions at <https://github.com/stevenliaotw/sltools>
pkg <- c("foreign", "Zelig", "reshape2",
         "ggplot2", "Amelia", "parallel",
         "plyr", "dplyr", "stringr", "sltools", 
         "tidyr", "gridExtra", "grid",
         "rworldmap", "rworldxtra", "RColorBrewer")
lapply(pkg, require, character.only = TRUE)

# set directory
MAIN_DIR <- "/Users/stevenliao/Dropbox/Research/No Reservations/replication"
# MAIN_DIR <- "replace above with your own directory here"

# create 2 folders named "data" and "out" within your main directory
# place replication data in the "data" folder
# the following code will load data from the "data" folder and output figures in the "out" folder
FIG_DIR <- paste(MAIN_DIR, "fig", sep = "/")
DATA_DIR <- paste(MAIN_DIR, "data", sep = "/")

# load data
load(paste(DATA_DIR, "merge.RData", sep = "/"))
load(paste(DATA_DIR, "reserve_mi.RData", sep = "/"))
load(paste(DATA_DIR, "reserve_mi_cs.RData", sep = "/"))

data.settlement <- read.csv(paste(DATA_DIR, "yuan.settlement.csv", sep = "/"), as.is = TRUE)
qfii <- read.csv(paste(DATA_DIR, "qfii_data.csv", sep = "/"), as.is = TRUE)
treasury <- read.csv(paste(DATA_DIR, "fredgraph.csv", sep = "/"), as.is = TRUE)
dollar.index <- read.csv(paste(DATA_DIR, "FRB_H10.csv", sep = "/"), as.is = TRUE)

# fix coding error for Japanese BSAs in 2012 and 2013
for(i in 1:length(data.mi$imputations)){
  data.mi$imputations[[i]][which(data.mi$imputations[[i]]$iso3c == "JPN" & data.mi$imputations[[i]]$year >= 2012),]$swap <- 1
}

for(i in 1:length(data.mi.cs)){
  data.mi.cs[[i]][which(data.mi.cs[[i]]$iso3c == "JPN"),]$swap <- 1
}

# set parameters
set.seed(1234)
mi.index <- c(1:10)

# create dataframe excluding SARs and TWN
nonUN.vec <- c("MAC", "HKG", "TWN")
data.mi.UN <- llply(data.mi$imputations, function(x) x[!(x$iso3c %in% nonUN.vec),])
data.mi.UN.cs <- llply(data.mi.cs, function(x) x[!(x$iso3c %in% nonUN.vec),])

# create dataframe excluding the US
nonUS.vec <- c("USA")
data.mi.nonUS <- llply(data.mi$imputations, function(x) x[!(x$iso3c %in% nonUS.vec),])

# set year as factor
for(i in 1:10){
  data.mi.UN[[i]]$year <- as.factor(data.mi.UN[[i]]$year)
}


###############################################################
## Fig 1. Countries Holding RMB Reserves in 2014. 
###############################################################
# subset
data.2013 <- filter(data, year == 2013)

# rworldmap
map.2013 <- joinCountryData2Map(data.2013, joinCode = "ISO3", nameJoinColumn = "iso3c", verbose = TRUE)

# create bins
#breaks <- c(0, 0.001, 0.01, 0.1, 1, 10, 45, 80)

# plot choropleth with rmb lead variable
mapDevice(device = "pdf", width = 7, file = paste(FIG_DIR, "chor_rmbres_2014.pdf", sep = "/"))
mapCountryData(map.2013, 
               nameColumnToPlot = "rmbres_lead", # thus plotting 2014 rmbres 
               mapRegion = "world",
               numCats = 2,
               catMethod = "categorical" ,
               missingCountryCol = "white",
               mapTitle = "", 
               addLegend = FALSE, 
               colourPalette = c("white", "darkred"),
               lwd = 1)
dev.off()


###############################################################
## Fig 2. QFII Expansion, 2003-2014. 
###############################################################
# create variable
qfii <- qfii %>%
  mutate(quota = quota/1000) # unit billion

# plot
p.qfii.1 <- ggplot(qfii, aes(x = year, y = quota)) + 
  geom_line(size = 0.6) + 
  xlab("") + ylab("Aggregate Quota (billion nominal USD)\n") + 
  scale_x_continuous(breaks = seq(2000, 2014, 1)) + 
  theme_bw() + theme(axis.title.y = element_text(size = 11),
                     plot.margin=unit(c(0.2, 0.2, -0.5, 0.2), "cm")
                     #axis.ticks = element_blank(), axis.text.x = element_blank()
                     ) 
p.qfii.2 <- ggplot(qfii, aes(x = year, y = institutions)) + 
  geom_line(size = 0.6) + 
  xlab("\nYear") + ylab("Number of Institutions (total QFIIs)\n") + 
  scale_x_continuous(breaks = seq(2000, 2014, 1)) + 
  scale_y_continuous(breaks = seq(0, 300, by = 50)) + 
  theme_bw() + theme(axis.title.x = element_text(size = 12),
                     axis.title.y = element_text(size = 11),
                     plot.margin=unit(c(0.2, 0.2, 0.2, 0.2), "cm")) 

pdf(paste(FIG_DIR, "qfii.pdf", sep = "/"), height = 6, width = 8)
grid.arrange(p.qfii.1, p.qfii.2, ncol=1)
dev.off()


###############################################################
## Fig 3. RMB Settlement of Cross-Border Trade.
###############################################################
# create variable
data.settlement$trade_set_per1 <- data.settlement$trade_set_per*100

# set limits
ylim1 <- max(data.settlement$trade_set_tot)
ylim2 <- max(data.settlement$trade_set_per1)

# plot
p.ts1 <- ggplot(data.settlement, 
                aes(x = time, y = trade_set_tot)) + 
  geom_bar(stat = "identity", width = .5) + 
  xlab("") + ylab("RMB-based Trade (USD billions)\n") + 
  scale_y_continuous(limits = c(0, ylim1)#, 
                     #breaks=c(0,25,50,75,100,125)
                     ) + 
  theme_bw() + theme(axis.title.y = element_text(size = 12),
                     axis.text.x  = element_text(angle = 30, 
                                                 vjust = 0.5, 
                                                 size = 11),
                     plot.margin = unit(c(0.2, 0.2, -0.2, 0.2), "cm")) 
p.ts2 <- ggplot(data.settlement, aes(x = time, y = trade_set_per1)) + 
  geom_bar(stat = "identity", width = .5) + 
  xlab("") + ylab("% of China's Total Trade\n") + 
  scale_y_continuous(limits = c(0, ylim2)#, 
                     #breaks=c(0,2,4,6,8,10,12)
                     ) + 
  theme_bw() + theme(axis.title.y = element_text(size = 12),
                     axis.text.x  = element_text(angle = 30, 
                                                 vjust = 0.5, 
                                                 size = 11),
                     plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), "cm")) 
pdf(paste(FIG_DIR, "trade_settle.pdf", sep = "/"), height = 6, width = 10)
grid.arrange(p.ts1, p.ts2, ncol=1)
dev.off()

###############################################################
## Fig 4. 5-Year Treasury and the Dollar’s Exchange Rate, 1973-2014 (monthly).
###############################################################
# load packages
library(lubridate)
library(scales)

# clean data for 5-year treasury bill
treas.df <- treasury[-c(1:10),]
names(treas.df) <- c("date", "treasury")
treas.df$date <- as.Date(parse_date_time(treas.df$date, "ymd"))
treas.df$treasury <- as.numeric(treas.df$treasury)
treas.df <- filter(treas.df, date >= as.Date("1973-01-01") & date <= as.Date("2014-12-01"))

# clean data for Dollar index
dollar.df <- dollar.index[-c(1:5),]
names(dollar.df) <- c("date", "dollar_index")
dollar.df$date <- as.Date(parse_date_time(dollar.df$date, "ym"))
dollar.df$dollar_index <- as.numeric(dollar.df$dollar_index)
dollar.df <- filter(dollar.df, date >= as.Date("1973-01-01") & date <= as.Date("2014-12-01"))

# plot
p.treas <- ggplot(treas.df, aes(x = date, y = treasury)) + 
  geom_line(size = 0.3) + 
  xlab("") + ylab("5-year Treasury Constant Monthly Rate\n") + 
  scale_x_date(breaks = date_breaks("2 years"),
               labels = date_format("%Y"),
               limits = c(as.Date("1974-01-01"), as.Date("2014-12-01"))) + 
  theme_bw() + theme(axis.title.y = element_text(size = 10),
                     axis.text.x  = element_text(angle = 60, 
                                                 vjust = 0.5, 
                                                 size = 10),
                     plot.margin=unit(c(0.3, 0.3, -0.5, 0.3), "cm")
                     #axis.ticks = element_blank(), axis.text.x = element_blank()
                     ) 

p.dollar <- ggplot(dollar.df, aes(x = date, y = dollar_index)) + 
  geom_line(size = 0.3) + 
  scale_x_date(breaks = date_breaks("2 years"),
               labels = date_format("%Y"),
               limits = c(as.Date("1974-01-01"), as.Date("2014-12-01"))) + 
  xlab("\nYear") + ylab("Nominal Major Currencies Dollar Index\n") + 
  theme_bw() + theme(axis.title.y = element_text(size = 10),
                     axis.text.x  = element_text(angle = 60, 
                                                 vjust = 0.5, 
                                                 size = 10),
                     plot.margin=unit(c(0.3, 0.3, 0.3, 0.3), "cm")) 

pdf(paste(FIG_DIR, "treas_dollar.pdf", sep = "/"), height = 6, width = 8)
grid.arrange(p.treas, p.dollar, ncol=1)
dev.off()


###############################################################
## Setup model specifications
###############################################################
# create vectors of control covariates
control.covars <- c("swap", "polity2", "log(gdp_2005_wdi)", "log(gdp_pc_2005_wdi)", "dist_k", "irto_cumsum_cbi",
                    "as.factor(year)")
print(control.covars)

## IO models
# create vectors of covariates
io.covars <- c("ipd_us",
               "ipd_cn",
               "s3un_us", 
               "s3un_cn")
               
# create lists of formulas
io.formula <- llply(io.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars), collapse = " + "))))
names(io.formula) <- io.covars
print(io.formula)

# for cross-sectional dropping year-fixed effects
io.formula.cs <- llply(io.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars[!grepl("year", control.covars)]), collapse = " + "))))
names(io.formula.cs) <- io.covars
print(io.formula.cs)


## TN Models (previously names TI model)
# create vectors of model covariates
ti.covars <- c("log(pcnimpdep)", "log(pcnimpdepgdp)")

# create lists of formulas
ti.formula <- llply(ti.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars), collapse = " + "))))
names(ti.formula) <- ti.covars
print(ti.formula)

# for cross-sectional dropping year-fixed effects
ti.formula.cs <- llply(ti.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars[!grepl("year", control.covars)]), collapse = " + "))))
names(ti.formula.cs) <- ti.covars
print(ti.formula.cs)


## OP models (previously named OC model)
# create vectors of covariates
oc.covars <- c("log(res_imp)",
                  #"log(res_debt)",
                  "log(res_gdp_cur)"#, 
                  #"res_gdp_change", 
                  #"res_gdp_change_08_per"
                  )
# create lists of formulas
oc.formula <- llply(oc.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars), collapse = " + "))))
names(oc.formula) <- oc.covars
print(oc.formula)

# for cross-sectional dropping year-fixed effects
oc.formula.cs <- llply(oc.covars, function(x) 
  formula(paste("rmbres_lead ~ ", paste(c(x, control.covars[!grepl("year", control.covars)]), collapse = " + "))))
names(oc.formula.cs) <- oc.covars
print(oc.formula.cs)

## IC models (Previously named QPQ model)
# create vectors of covariates
qpq.covars <- c("log(pcnexpdep)", "pfdidep", "log(pcnexpdepgdp)")

# create lists of formulas
qpq.formula <- list()
qpq.formula$`log(pcnexpdep)` <- formula(paste("rmbres_lead ~ ", 
                                  paste(c(qpq.covars[1:2], control.covars), 
                                        collapse = " + ")))
qpq.formula$`log(pcnexpdepgdp)` <- formula(paste("rmbres_lead ~ ", 
                                    paste(c(qpq.covars[c(3, 2)], control.covars), 
                                          collapse = " + ")))
print(qpq.formula)

# for cross-sectional dropping year-fixed effects
qpq.formula.cs <- list()
qpq.formula.cs$`log(pcnexpdep)` <- formula(paste("rmbres_lead ~ ", 
                                     paste(c(qpq.covars[1:2], 
                                             control.covars[!grepl("year", control.covars)]), 
                                     collapse = " + ")))
qpq.formula.cs$`log(pcnexpdepgdp)` <- formula(paste("rmbres_lead ~ ", 
                                       paste(c(qpq.covars[c(3, 2)], 
                                               control.covars[!grepl("year", control.covars)]), 
                                       collapse = " + ")))
print(qpq.formula.cs)

## Full models
# create vectors of covariates
full.covars.us <- c("log(pcnimpdep)", "log(res_imp)", "log(pcnexpdep)", "pfdidep", "ipd_us")
full.covars.cn <- c("log(pcnimpdep)", "log(res_imp)", "log(pcnexpdep)", "pfdidep", "ipd_cn")

# create vectors of full model covariates
full.formula.us <- formula(paste("rmbres_lead ~ ", paste(c(full.covars.us, control.covars), collapse = " + ")))
print(full.formula.us)

full.formula.cn <- formula(paste("rmbres_lead ~ ", paste(c(full.covars.cn, control.covars), collapse = " + ")))
print(full.formula.cn)

full.formula.us.cs <- formula(paste("rmbres_lead ~ ", paste(c(full.covars.us, control.covars[!grepl("year", control.covars)]), collapse = " + ")))
print(full.formula.us.cs)

full.formula.cn.cs <- formula(paste("rmbres_lead ~ ", paste(c(full.covars.cn, control.covars[!grepl("year", control.covars)]), collapse = " + ")))
print(full.formula.cn.cs)


###############################################################
## Fit logit models
###############################################################
## IO (UN countries only: exclude MAC, HKG, TWN)
# using ideal point distance for USA
io.logit.us <- zelig(io.formula$ipd_us,
                   model = "logit", 
                   data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                             data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                   robust = TRUE)
io.logit.us.sum <- summary(io.logit.us)
io.logit.us.sum

# using ideal point distance for CHN
io.logit.cn <- zelig(io.formula$ipd_cn,
                   model = "logit", 
                   data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                             data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                   robust = TRUE)
io.logit.cn.sum <- summary(io.logit.cn)
io.logit.cn.sum

## TI
ti.logit <- zelig(ti.formula$`log(pcnimpdep)`,
                  model = "logit", 
                  data = mi(data.mi$imputations[[1]], data.mi$imputations[[2]], data.mi$imputations[[3]], data.mi$imputations[[4]], data.mi$imputations[[5]],
                            data.mi$imputations[[6]], data.mi$imputations[[7]], data.mi$imputations[[8]], data.mi$imputations[[9]], data.mi$imputations[[10]]), 
                  robust = TRUE)
ti.logit.sum <- summary(ti.logit)
ti.logit.sum

## OC
oc.logit <- zelig(oc.formula$`log(res_imp)`,
                  model = "logit",
                  data = mi(data.mi$imputations[[1]], data.mi$imputations[[2]], data.mi$imputations[[3]], data.mi$imputations[[4]], data.mi$imputations[[5]],
                            data.mi$imputations[[6]], data.mi$imputations[[7]], data.mi$imputations[[8]], data.mi$imputations[[9]], data.mi$imputations[[10]]), 
                  robust = TRUE)
oc.logit.sum <- summary(oc.logit)
oc.logit.sum

## QPQ
qpq.logit <- zelig(qpq.formula$`log(pcnexpdep)`,
                   model = "logit", 
                   data = mi(data.mi$imputations[[1]], data.mi$imputations[[2]], data.mi$imputations[[3]], data.mi$imputations[[4]], data.mi$imputations[[5]],
                             data.mi$imputations[[6]], data.mi$imputations[[7]], data.mi$imputations[[8]], data.mi$imputations[[9]], data.mi$imputations[[10]]), 
                   robust = TRUE)
qpq.logit.sum <- summary(qpq.logit)
qpq.logit.sum


## Full model (UN countries only: exclude MAC, HKG, TWN)
# using ideal point distance for USA
full.logit.us <- zelig(full.formula.us,
                      model = "logit", 
                      data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                      robust = TRUE)
full.logit.us.sum <- summary(full.logit.us)
full.logit.us.sum

# using ideal point distance for CHN
full.logit.cn <- zelig(full.formula.cn,
                      model = "logit", 
                      data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                      robust = TRUE)
full.logit.cn.sum <- summary(full.logit.cn)
full.logit.cn.sum


###############################################################
## Fig 5. Results from Fitted Logistic Regression Models. 
###############################################################
# make list of fitted model results
#alpha <- 0.05
#Multiplier <- qnorm(1 - alpha / 2)
models <- list(io.logit.us, io.logit.cn, ti.logit, oc.logit, qpq.logit, full.logit.us, full.logit.cn)
model.names <- c("IO (US)", "IO (CN)", "TN", "OP", "IC", "Full (US)", "Full (CN)")
var.order <- c("ipd_us", "ipd_cn", "log(pcnimpdep)", "log(res_imp)", "log(pcnexpdep)", "pfdidep", "swap", "polity2", "log(gdp_2005_wdi)", "log(gdp_pc_2005_wdi)", "dist_k", "irto_cumsum_cbi")
var.names <- c("UNGA IPD (US)",
               "UNGA IPD (CN)",
               "Log (Partner Import Dep., total imports)",
               "Log (Reserves in Months of Imports)",
               "Log (Partner Export Dep., total exports)",
               "Partner FDI Inflow Dep.",
               "Bilateral Swap Agreement",
               "Regime Type",
               "Log (GDP)",
               "Log (GDP per Capita)",
               "Distance (thousand km)",
               "Cumulative Irregular CB Turnovers")

coef.plot <- coefficientPlot(models, alpha = 0.05, 
                             model.names = model.names, 
                             var.order = var.order, 
                             var.names = var.names, 
                             omit.intercept = TRUE, omit.factor = TRUE)

pdf(file = paste(FIG_DIR, "coefplot.pdf", sep = "/"), width = 12, height = 8)
coef.plot
dev.off()



###############################################################
## FIG 6. Substantive Effects
###############################################################
## ipd_us
# extract parameters
coef.ipd.us <- ldply(full.logit.us, function(x) coef(x))
vcov.ipd.us <- llply(full.logit.us, function(x) vcov(x))

# Set x values  
ipd.us.range <- with(data, seq(min(ipd_us, na.rm = TRUE), max(ipd_us, na.rm = TRUE), length.out = 50))
typical.ipd.us.cty.year <- setx(full.logit.us, ipd_us = ipd.us.range, year = 2012)

# simultate predicted probabilities
ipd.us.qi <- computeLogitEV(data = data.mi.UN, x = typical.ipd.us.cty.year, coef = coef.ipd.us, vcov = vcov.ipd.us, num = 100)

# check dimensions
dim(ipd.us.qi$ev)

# merge with set values
ipd.us.ev <- as.data.frame(ipd.us.qi$ev)
ipd.us.ev <- cbind(ipd.us.ev, ipd.us.range)

# reshpae
ipd.us.ev.long <- ipd.us.ev %>% 
  gather(draw, value, -ipd.us.range) %>%
  arrange(ipd.us.range)


## ipd_cn
# extract parameters
coef.ipd.cn <- ldply(full.logit.cn, function(x) coef(x))
vcov.ipd.cn <- llply(full.logit.cn, function(x) vcov(x))

# Set x values  
ipd.cn.range <- with(data, seq(min(ipd_cn, na.rm = TRUE), max(ipd_cn, na.rm = TRUE), length.out = 50))
typical.ipd.cn.cty.year <- setx(full.logit.cn, ipd_cn = ipd.cn.range, year = 2012)

# simultate predicted probabilities
ipd.cn.qi <- computeLogitEV(data = data.mi.UN, x = typical.ipd.cn.cty.year, coef = coef.ipd.cn, vcov = vcov.ipd.cn, num = 100)

# check dimensions
dim(ipd.cn.qi$ev)

# merge with set values
ipd.cn.ev <- as.data.frame(ipd.cn.qi$ev)
ipd.cn.ev <- cbind(ipd.cn.ev, ipd.cn.range)

# reshpae
ipd.cn.ev.long <- ipd.cn.ev %>% 
  gather(draw, value, -ipd.cn.range) %>%
  arrange(ipd.cn.range)

# ribbon plots
# ipd_us
# calculate mean, upper/lower bounds
s.ipd.us.df <- data.frame(x = ipd.us.range,
                          mean = apply(ipd.us.ev, 1, mean),
                          lwr = apply(ipd.us.ev, 1, quantile, probs = 0.025),
                          upr = apply(ipd.us.ev, 1, quantile, probs = 0.975))

# plot
p.ipd.us.ribbon <- ggplot(s.ipd.us.df, 
                          aes(x = x, y = mean)) + 
  geom_ribbon(aes(ymin = lwr,
                  ymax = upr), 
              alpha = 0.2) + 
  geom_line(size = 1) + 
  scale_x_continuous("\nUN Ideal Point Distance with USA", expand = c(0, 0)) + 
  scale_y_continuous("Predicted Probability of RMB Reserves\n", 
                     limits = c(0, 1), expand = c(0, 0)) +
  theme_bw() + geom_vline(xintercept = c(mean(data$ipd_us, na.rm = TRUE) - sd(data$ipd_us, na.rm = TRUE),
                                         mean(data$ipd_us, na.rm = TRUE) + sd(data$ipd_us, na.rm = TRUE)),
                          colour = I(hsv(0/12, 7/12, 7/12)), linetype = "longdash")

# ipd_cn
# calculate mean, upper/lower bounds
s.ipd.cn.df <- data.frame(x = ipd.cn.range,
                          mean = apply(ipd.cn.ev, 1, mean),
                          lwr = apply(ipd.cn.ev, 1, quantile, probs = 0.025),
                          upr = apply(ipd.cn.ev, 1, quantile, probs = 0.975))

# plot
p.ipd.cn.ribbon <- ggplot(s.ipd.cn.df, 
                          aes(x = x, y = mean)) + 
  geom_ribbon(aes(ymin = lwr,
                  ymax = upr), 
              alpha = 0.2) + 
  geom_line(size = 1) + 
  scale_x_continuous("\nUN Ideal Point Distance with China", expand = c(0, 0)) + 
  scale_y_continuous("Predicted Probability of RMB Reserves\n", 
                     limits = c(0, 1), expand = c(0, 0)) +
  theme_bw() + geom_vline(xintercept = c(mean(data$ipd_cn, na.rm = TRUE) - sd(data$ipd_cn, na.rm = TRUE),
                                         mean(data$ipd_cn, na.rm = TRUE) + sd(data$ipd_cn, na.rm = TRUE)),
                          colour = I(hsv(0/12, 7/12, 7/12)), linetype = "longdash")

pdf(paste(FIG_DIR, "effects_ribbon.pdf", sep = "/"), width = 10, height = 4)  
grid.arrange(p.ipd.us.ribbon, p.ipd.cn.ribbon,
             ncol = 2)
dev.off()


###############################################################
## Fit robustness models: alternative measures
###############################################################
## IO (UN countries only: exclude MAC, HKG, TWN)
# USA
io.logit.robust.s3un.us <- zelig(rmbres_lead ~ log(pcnimpdep) + log(res_imp) + log(pcnexpdep) + 
                                 pfdidep + s3un_us + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                                 dist_k + irto_cumsum_cbi + as.factor(year),
                                 model = "logit", 
                                 data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                           data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                                 robust = TRUE)
io.logit.robust.s3un.us.sum <- summary(io.logit.robust.s3un.us)
io.logit.robust.s3un.us.sum

# CHN
io.logit.robust.s3un.cn <- zelig(rmbres_lead ~ log(pcnimpdep) + log(res_imp) + log(pcnexpdep) + 
                                 pfdidep + s3un_cn + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                                 dist_k + irto_cumsum_cbi + as.factor(year),
                                 model = "logit", 
                                 data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                           data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                                 robust = TRUE)
io.logit.robust.s3un.cn.sum <- summary(io.logit.robust.s3un.cn)
io.logit.robust.s3un.cn.sum

## TI
ti.logit.robust <- zelig(rmbres_lead ~ log(pcnimpdepgdp) + log(res_imp) + log(pcnexpdep) + 
                         pfdidep + ipd_us + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                         dist_k + irto_cumsum_cbi + as.factor(year),
                         model = "logit", 
                         data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                   data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                         robust = TRUE)
ti.logit.robust.sum <- summary(ti.logit.robust)
ti.logit.robust.sum

## OC
oc.logit.robust <- zelig(rmbres_lead ~ log(pcnimpdep) + log(res_gdp_cur) + log(pcnexpdep) + 
                         pfdidep + ipd_us + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                         dist_k + irto_cumsum_cbi + as.factor(year),
                         model = "logit", 
                         data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                   data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                         robust = TRUE)
oc.logit.robust.sum <- summary(oc.logit.robust)
oc.logit.robust.sum

## QPQ
qpq.logit.robust <- zelig(rmbres_lead ~ log(pcnimpdep) + log(res_imp) + log(pcnexpdepgdp) + 
                            pfdidep + ipd_us + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                            dist_k + irto_cumsum_cbi + as.factor(year),
                          model = "logit", 
                          data = mi(data.mi.UN[[1]], data.mi.UN[[2]], data.mi.UN[[3]], data.mi.UN[[4]], data.mi.UN[[5]],
                                    data.mi.UN[[6]], data.mi.UN[[7]], data.mi.UN[[8]], data.mi.UN[[9]], data.mi.UN[[10]]), 
                          robust = TRUE)
qpq.logit.robust.sum <- summary(qpq.logit.robust)
qpq.logit.robust.sum

###############################################################
## Fit robustness models: cross-sectional
###############################################################
# collapse into cross-sectional analysis
# idp with USA
io.logit.fm.us.cs <- zelig(io.formula.cs$ipd_us,
                           model = "logit", 
                           data = mi(data.mi.UN.cs[[1]], data.mi.UN.cs[[2]], 
                                     data.mi.UN.cs[[3]], data.mi.UN.cs[[4]], 
                                     data.mi.UN.cs[[5]], data.mi.UN.cs[[6]], 
                                     data.mi.UN.cs[[7]], data.mi.UN.cs[[8]], 
                                     data.mi.UN.cs[[9]], data.mi.UN.cs[[10]]), 
                           robust = TRUE)
io.logit.fm.us.cs.sum <- summary(io.logit.fm.us.cs)
io.logit.fm.us.cs.sum

# idp with CHN
io.logit.fm.cn.cs <- zelig(io.formula.cs$ipd_cn,
                           model = "logit", 
                           data = mi(data.mi.UN.cs[[1]], data.mi.UN.cs[[2]], 
                                     data.mi.UN.cs[[3]], data.mi.UN.cs[[4]], 
                                     data.mi.UN.cs[[5]], data.mi.UN.cs[[6]], 
                                     data.mi.UN.cs[[7]], data.mi.UN.cs[[8]], 
                                     data.mi.UN.cs[[9]], data.mi.UN.cs[[10]]), 
                           robust = TRUE)
io.logit.fm.cn.cs.sum <- summary(io.logit.fm.cn.cs)
io.logit.fm.cn.cs.sum

# using pcnimpdep
ti.logit.fm.cs <- zelig(ti.formula.cs$`log(pcnimpdep)`,
                        model = "logit", 
                        data = mi(data.mi.cs[[1]], data.mi.cs[[2]], 
                                  data.mi.cs[[3]], data.mi.cs[[4]], 
                                  data.mi.cs[[5]], data.mi.cs[[6]], 
                                  data.mi.cs[[7]], data.mi.cs[[8]], 
                                  data.mi.cs[[9]], data.mi.cs[[10]]), 
                        robust = TRUE)
ti.logit.fm.cs.sum <- summary(ti.logit.fm.cs)
ti.logit.fm.cs.sum

# using res_imp
oc.logit.fm.cs <- zelig(oc.formula.cs$`log(res_imp)`,
                        model = "logit", 
                        data = mi(data.mi.cs[[1]], data.mi.cs[[2]], 
                                  data.mi.cs[[3]], data.mi.cs[[4]], 
                                  data.mi.cs[[5]], data.mi.cs[[6]], 
                                  data.mi.cs[[7]], data.mi.cs[[8]], 
                                  data.mi.cs[[9]], data.mi.cs[[10]]), 
                        robust = TRUE)
oc.logit.fm.cs.sum <- summary(oc.logit.fm.cs)
oc.logit.fm.cs.sum

# using log(pcnexpdep)
qpq.logit.fm.cs <- zelig(qpq.formula.cs$`log(pcnexpdep)`,
                         model = "logit", 
                         data = mi(data.mi.cs[[1]], data.mi.cs[[2]], 
                                   data.mi.cs[[3]], data.mi.cs[[4]], 
                                   data.mi.cs[[5]], data.mi.cs[[6]], 
                                   data.mi.cs[[7]], data.mi.cs[[8]], 
                                   data.mi.cs[[9]], data.mi.cs[[10]]), 
                         robust = TRUE)
qpq.logit.fm.cs.sum <- summary(qpq.logit.fm.cs)
qpq.logit.fm.cs.sum

## full model in cross-sectional 
# USA
full.logit.fm.us.cs <- zelig(full.formula.us.cs,
                             model = "logit", 
                             data = mi(data.mi.UN.cs[[1]], data.mi.UN.cs[[2]], 
                                       data.mi.UN.cs[[3]], data.mi.UN.cs[[4]], 
                                       data.mi.UN.cs[[5]], data.mi.UN.cs[[6]], 
                                       data.mi.UN.cs[[7]], data.mi.UN.cs[[8]], 
                                       data.mi.UN.cs[[9]], data.mi.UN.cs[[10]]), 
                             robust = TRUE)
full.logit.fm.us.cs.sum <- summary(full.logit.fm.us.cs)
full.logit.fm.us.cs.sum

# CHN
full.logit.fm.cn.cs <- zelig(full.formula.cn.cs,
                             model = "logit", 
                             data = mi(data.mi.UN.cs[[1]], data.mi.UN.cs[[2]], 
                                       data.mi.UN.cs[[3]], data.mi.UN.cs[[4]], 
                                       data.mi.UN.cs[[5]], data.mi.UN.cs[[6]], 
                                       data.mi.UN.cs[[7]], data.mi.UN.cs[[8]], 
                                       data.mi.UN.cs[[9]], data.mi.UN.cs[[10]]), 
                             robust = TRUE)
full.logit.fm.cn.cs.sum <- summary(full.logit.fm.cn.cs)
full.logit.fm.cn.cs.sum

###############################################################
## Fit robustness models: conditional logit
###############################################################
library(survival)

## IO
# USA
full.clogit.us <- llply(data.mi.UN, function(x) {
  return(clogit(rmbres_lead ~ log(pcnimpdep) + log(res_imp) + log(pcnexpdep) + 
                  pfdidep + ipd_us + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                  dist_k + irto_cumsum_cbi + strata(year),  data = x))
})
full.clogit.us.summary <- combineCLogitMI(full.clogit.us)

# CHN
full.clogit.cn <- llply(data.mi.UN, function(x) {
  return(clogit(rmbres_lead ~ log(pcnimpdep) + log(res_imp) + log(pcnexpdep) + 
                  pfdidep + ipd_cn + swap + polity2 + log(gdp_2005_wdi) + log(gdp_pc_2005_wdi) + 
                  dist_k + irto_cumsum_cbi + strata(year),  data = x))
})
full.clogit.cn.summary <- combineCLogitMI(full.clogit.cn)




###############################################################
## Extract Results for Tables
###############################################################
# IO
io.logit.us.texreg <- texregMI(names = rownames(coef(io.logit.us.sum)), 
                               coef = coef(io.logit.us.sum)[,1], 
                               se = coef(io.logit.us.sum)[,2], 
                               pval = coef(io.logit.us.sum)[,4], 
                               n = length(io.logit.us[[1]]$data$iso3c),
                               loglik = mean(sapply(mi.index, function(x) io.logit.us[[x]]$deviance/-2)),
                               aic = mean(sapply(mi.index, function(x) io.logit.us[[x]]$aic)))
io.logit.cn.texreg <- texregMI(names = rownames(coef(io.logit.cn.sum)), 
                               coef = coef(io.logit.cn.sum)[,1], 
                               se = coef(io.logit.cn.sum)[,2], 
                               pval = coef(io.logit.cn.sum)[,4], 
                               n = length(io.logit.cn[[1]]$data$iso3c),
                               loglik = mean(sapply(mi.index, function(x) io.logit.cn[[x]]$deviance/-2)),
                               aic = mean(sapply(mi.index, function(x) io.logit.cn[[x]]$aic)))
io.logit.robust.s3un.us.texreg <- texregMI(names = rownames(coef(io.logit.robust.s3un.us.sum)), 
                                           coef = coef(io.logit.robust.s3un.us.sum)[,1], 
                                           se = coef(io.logit.robust.s3un.us.sum)[,2], 
                                           pval = coef(io.logit.robust.s3un.us.sum)[,4], 
                                           n = length(io.logit.robust.s3un.us[[1]]$data$iso3c),
                                           loglik = mean(sapply(mi.index, function(x) io.logit.robust.s3un.us[[x]]$deviance/-2)),
                                           aic = mean(sapply(mi.index, function(x) io.logit.robust.s3un.us[[x]]$aic)))
io.logit.robust.s3un.cn.texreg <- texregMI(names = rownames(coef(io.logit.robust.s3un.cn.sum)), 
                                           coef = coef(io.logit.robust.s3un.cn.sum)[,1], 
                                           se = coef(io.logit.robust.s3un.cn.sum)[,2], 
                                           pval = coef(io.logit.robust.s3un.cn.sum)[,4], 
                                           n = length(io.logit.robust.s3un.cn[[1]]$data$iso3c),
                                           loglik = mean(sapply(mi.index, function(x) io.logit.robust.s3un.cn[[x]]$deviance/-2)),
                                           aic = mean(sapply(mi.index, function(x) io.logit.robust.s3un.cn[[x]]$aic)))
io.logit.us.cs.texreg <- texregMI(names = rownames(coef(io.logit.fm.us.cs.sum)), 
                                  coef = coef(io.logit.fm.us.cs.sum)[,1], 
                                  se = coef(io.logit.fm.us.cs.sum)[,2], 
                                  pval = coef(io.logit.fm.us.cs.sum)[,4], 
                                  n = length(io.logit.fm.us.cs[[1]]$data$iso3c),
                                  loglik = mean(sapply(mi.index, function(x) io.logit.fm.us.cs[[x]]$deviance/-2)),
                                  aic = mean(sapply(mi.index, function(x) io.logit.fm.us.cs[[x]]$aic)))
io.logit.cn.cs.texreg <- texregMI(names = rownames(coef(io.logit.fm.cn.cs.sum)), 
                                  coef = coef(io.logit.fm.cn.cs.sum)[,1], 
                                  se = coef(io.logit.fm.cn.cs.sum)[,2], 
                                  pval = coef(io.logit.fm.cn.cs.sum)[,4], 
                                  n = length(io.logit.fm.cn.cs[[1]]$data$iso3c),
                                  loglik = mean(sapply(mi.index, function(x) io.logit.fm.cn.cs[[x]]$deviance/-2)),
                                  aic = mean(sapply(mi.index, function(x) io.logit.fm.cn.cs[[x]]$aic)))

# TI
ti.logit.texreg <- texregMI(names = rownames(coef(ti.logit.sum)), 
                            coef = coef(ti.logit.sum)[,1], 
                            se = coef(ti.logit.sum)[,2], 
                            pval = coef(ti.logit.sum)[,4], 
                            n = length(ti.logit[[1]]$data$iso3c),
                            loglik = mean(sapply(mi.index, function(x) ti.logit[[x]]$deviance/-2)),
                            aic = mean(sapply(mi.index, function(x) ti.logit[[x]]$aic)))
ti.logit.robust.texreg <- texregMI(names = rownames(coef(ti.logit.robust.sum)), 
                                   coef = coef(ti.logit.robust.sum)[,1], 
                                   se = coef(ti.logit.robust.sum)[,2], 
                                   pval = coef(ti.logit.robust.sum)[,4], 
                                   n = length(ti.logit.robust[[1]]$data$iso3c),
                                   loglik = mean(sapply(mi.index, function(x) ti.logit.robust[[x]]$deviance/-2)),
                                   aic = mean(sapply(mi.index, function(x) ti.logit.robust[[x]]$aic)))
ti.logit.cs.texreg <- texregMI(names = rownames(coef(ti.logit.fm.cs.sum)), 
                               coef = coef(ti.logit.fm.cs.sum)[,1], 
                               se = coef(ti.logit.fm.cs.sum)[,2], 
                               pval = coef(ti.logit.fm.cs.sum)[,4], 
                               n = length(ti.logit.fm.cs[[1]]$data$iso3c),
                               loglik = mean(sapply(mi.index, function(x) ti.logit.fm.cs[[x]]$deviance/-2)),
                               aic = mean(sapply(mi.index, function(x) ti.logit.fm.cs[[x]]$aic)))
# OC
oc.logit.texreg <- texregMI(names = rownames(coef(oc.logit.sum)), 
                            coef = coef(oc.logit.sum)[,1], 
                            se = coef(oc.logit.sum)[,2], 
                            pval = coef(oc.logit.sum)[,4], 
                            n = length(oc.logit[[1]]$data$iso3c),
                            loglik = mean(sapply(mi.index, function(x) oc.logit[[x]]$deviance/-2)),
                            aic = mean(sapply(mi.index, function(x) oc.logit[[x]]$aic)))
oc.logit.robust.texreg <- texregMI(names = rownames(coef(oc.logit.robust.sum)), 
                                   coef = coef(oc.logit.robust.sum)[,1], 
                                   se = coef(oc.logit.robust.sum)[,2], 
                                   pval = coef(oc.logit.robust.sum)[,4], 
                                   n = length(oc.logit.robust[[1]]$data$iso3c),
                                   loglik = mean(sapply(mi.index, function(x) oc.logit.robust[[x]]$deviance/-2)),
                                   aic = mean(sapply(mi.index, function(x) oc.logit.robust[[x]]$aic)))
oc.logit.cs.texreg <- texregMI(names = rownames(coef(oc.logit.fm.cs.sum)), 
                               coef = coef(oc.logit.fm.cs.sum)[,1], 
                               se = coef(oc.logit.fm.cs.sum)[,2], 
                               pval = coef(oc.logit.fm.cs.sum)[,4], 
                               n = length(oc.logit.fm.cs[[1]]$data$iso3c),
                               loglik = mean(sapply(mi.index, function(x) oc.logit.fm.cs[[x]]$deviance/-2)),
                               aic = mean(sapply(mi.index, function(x) oc.logit.fm.cs[[x]]$aic)))
# QPQ
qpq.logit.texreg <- texregMI(names = rownames(coef(qpq.logit.sum)), 
                             coef = coef(qpq.logit.sum)[,1], 
                             se = coef(qpq.logit.sum)[,2], 
                             pval = coef(qpq.logit.sum)[,4], 
                             n = length(qpq.logit[[1]]$data$iso3c),
                             loglik = mean(sapply(mi.index, function(x) qpq.logit[[x]]$deviance/-2)),
                             aic = mean(sapply(mi.index, function(x) qpq.logit[[x]]$aic)))
qpq.logit.robust.texreg <- texregMI(names = rownames(coef(qpq.logit.robust.sum)), 
                                    coef = coef(qpq.logit.robust.sum)[,1], 
                                    se = coef(qpq.logit.robust.sum)[,2], 
                                    pval = coef(qpq.logit.robust.sum)[,4], 
                                    n = length(qpq.logit.robust[[1]]$data$iso3c),
                                    loglik = mean(sapply(mi.index, function(x) qpq.logit.robust[[x]]$deviance/-2)),
                                    aic = mean(sapply(mi.index, function(x) qpq.logit.robust[[x]]$aic)))
qpq.logit.cs.texreg <- texregMI(names = rownames(coef(qpq.logit.fm.cs.sum)), 
                                coef = coef(qpq.logit.fm.cs.sum)[,1], 
                                se = coef(qpq.logit.fm.cs.sum)[,2], 
                                pval = coef(qpq.logit.fm.cs.sum)[,4], 
                                n = length(qpq.logit.fm.cs[[1]]$data$iso3c),
                                loglik = mean(sapply(mi.index, function(x) qpq.logit.fm.cs[[x]]$deviance/-2)),
                                aic = mean(sapply(mi.index, function(x) qpq.logit.fm.cs[[x]]$aic)))

# Full
full.logit.us.texreg <- texregMI(names = rownames(coef(full.logit.us.sum)), 
                                 coef = coef(full.logit.us.sum)[,1], 
                                 se = coef(full.logit.us.sum)[,2], 
                                 pval = coef(full.logit.us.sum)[,4], 
                                 n = length(full.logit.us[[1]]$data$iso3c),
                                 loglik = mean(sapply(mi.index, function(x) full.logit.us[[x]]$deviance/-2)),
                                 aic = mean(sapply(mi.index, function(x) full.logit.us[[x]]$aic)))
full.logit.cn.texreg <- texregMI(names = rownames(coef(full.logit.cn.sum)), 
                                 coef = coef(full.logit.cn.sum)[,1], 
                                 se = coef(full.logit.cn.sum)[,2], 
                                 pval = coef(full.logit.cn.sum)[,4], 
                                 n = length(full.logit.cn[[1]]$data$iso3c),
                                 loglik = mean(sapply(mi.index, function(x) full.logit.cn[[x]]$deviance/-2)),
                                 aic = mean(sapply(mi.index, function(x) full.logit.cn[[x]]$aic)))
full.logit.us.cs.texreg <- texregMI(names = rownames(coef(full.logit.fm.us.cs.sum)), 
                                    coef = coef(full.logit.fm.us.cs.sum)[,1], 
                                    se = coef(full.logit.fm.us.cs.sum)[,2], 
                                    pval = coef(full.logit.fm.us.cs.sum)[,4], 
                                    n = length(full.logit.fm.us.cs[[1]]$data$iso3c),
                                    loglik = mean(sapply(mi.index, function(x) full.logit.fm.us.cs[[x]]$deviance/-2)),
                                    aic = mean(sapply(mi.index, function(x) full.logit.fm.us.cs[[x]]$aic)))
full.logit.cn.cs.texreg <- texregMI(names = rownames(coef(full.logit.fm.cn.cs.sum)), 
                                    coef = coef(full.logit.fm.cn.cs.sum)[,1], 
                                    se = coef(full.logit.fm.cn.cs.sum)[,2], 
                                    pval = coef(full.logit.fm.cn.cs.sum)[,4], 
                                    n =  length(full.logit.fm.cn.cs[[1]]$data$iso3c),
                                    loglik = mean(sapply(mi.index, function(x) full.logit.fm.cn.cs[[x]]$deviance/-2)),
                                    aic = mean(sapply(mi.index, function(x) full.logit.fm.cn.cs[[x]]$aic)))

##############################################################################
## Figure B.1. Covariate Correlation Matrix 
##############################################################################
library(polycor)
library(reshape2)
library(ggplot2)
library(mgcv)

# subset vars
plot.df <- data.mi$imputations[[1]] %>%
  select(rmbres_lead,
         ipd_us,
         ipd_cn,
         s3un_us,
         s3un_cn,
         pcnimpdep, #ppgd_other, 
         pcnimpdepgdp,
         res_imp, 
         res_gdp_cur,
         pcnexpdep, 
         pcnexpdepgdp,
         pfdidep, 
         #scofull, 
         swap, 
         polity2, 
         gdp_2005_wdi, 
         gdp_pc_2005_wdi, 
         dist_k, 
         irto_cumsum_cbi) 
names(plot.df) <- c("RMB Reserves",
             "UNGA Ideal Point Distance (USA)",
             "UNGA Ideal Point Distance (China)",
             "UNGA Voting Affinity (USA)",
             "UNGA Voting Affinity (China)",
             "Partner Import Dependence (total imports)",
             "Partner Import Dependence (current GDP)",
             #"Debt in Other Currency",
             "Reserves in Months of Imports",
             "Reserves over Current GDP",
             "Partner Export Dependence (total exports)",
             "Partner Export Dependence (current GDP)",
             "Partner FDI Inflow Dependence",
             #"SCO",
             "BSA",
             "Regime Type",
             "GDP",
             "GDP per Capita",
             "Distance (thousand km)",
             "Cumulative Irregular CB Turnovers")

plot.df <- melt(hetcor(plot.df, use = "pairwise.complete.obs")$correlations)

plot.df <- plot.df %>%
  arrange(value)

# sort
plot.df$Var1 <- reorder(plot.df$Var1, plot.df$value)
plot.df$Var2 <- reorder(plot.df$Var2, plot.df$value)

# plot
p <- ggplot(data = plot.df, aes(x = Var1, y = Var2, fill = value))
p <- p + geom_tile()
p <- p + scale_fill_gradient2(name = "Correlation", breaks = seq(-1, 1, by = .25),
                              space = "Lab")
p <- p + guides(fill = guide_colorbar(barwidth = .75, ticks = FALSE))
p <- p + labs(x = NULL, y = NULL)
p <- p + theme(axis.text.x = element_text(angle = 60, hjust = 1))
pdf(file = paste(FIG_DIR, "corFig.pdf", sep = "/"), width = 8, height = 8)
p
dev.off()

##############################################################################
## TABLE C.1. Results from Main Logistic Regression Models Fitted to Panel Data
##############################################################################
texreg(list(io.logit.us.texreg, 
            io.logit.cn.texreg, 
            ti.logit.texreg, 
            oc.logit.texreg, 
            qpq.logit.texreg, 
            full.logit.us.texreg,
            full.logit.cn.texreg), 
       custom.model.names = c("IO (US)", "IO (CN)", 
                              "TN", "OP", "IC", 
                              "Full (US)", "Full (CN)"), 
       custom.coef.names = c("Intercept",  
                             "UNGA Ideal Point Distance with the USA",
                             "BSA",
                             "Regime Type",
                             "Log(GDP)",
                             "Log(GDP per Capita)",
                             "Distance (thousand km)",
                             "Cumulative Irregular CB Turnovers",
                             rep(NA, 4),
                             "UNGA Ideal Point Distance with the China",
                             "Log(Partner Import Dependence)",
                             "Log(Reserves in Months of Imports)",
                             "Log(Partner Export Dependence)",
                             "Partner FDI Inflow Dependence"),
       reorder.coef = c(1, 2, c(9:13), c(3:8)),
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol="\\wedge",
       digits = 3, omit.coef = "factor",
       fontsize = "footnotesize",
       single.row = FALSE, float.pos = "!ht", sideways = TRUE,
       caption = "Results for Main Logistic Regression Models Fitted to Panel Data. Dichotomous variable for RMB reserves in year $t +  1$ as the dependent variable. HAC standard errors are in parentheses. Year fixed effects included in model but results omitted to facilitate presentation.", 
       label = "tb:fm", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)


##############################################################################
## TABLE C.2. Results from Main Logistic Regression Models Fitted to Panel Data
##############################################################################
texreg(list(io.logit.robust.s3un.us.texreg, 
            io.logit.robust.s3un.cn.texreg,
            ti.logit.robust.texreg, 
            oc.logit.robust.texreg, 
            qpq.logit.robust.texreg,
            full.clogit.us.summary$texreg, 
            full.clogit.cn.summary$texreg), 
       custom.model.names = c("Full (Affinity US)", "Full (Affinity CN)",
                              "Full (alt. TN)", "Full (alt. OP)", "Full (alt. IC)", 
                              "Full (clogit US)", "Full (clogit CN)"),
       custom.coef.names = c("Intercept", 
                             "Log(Partner Import Dependence)",
                             "Log(Reserves in Months of Imports)",
                             "Log(Partner Export Dependence)",
                             "Partner FDI Inflow Dependence",
                             "UNGA Voting Affinity with the US",
                             "Bilateral Swap Agreements",
                             "Regime Type",
                             "Log(GDP)",
                             "Log(GDP per Capita)",
                             "Distance",
                             "Cumulative Irregular CB Turnovers",
                             rep(NA, 4),
                             "UNGA Voting Affinity with China",
                             "Log(Partner Import Dependence, GDP)",
                             "UNGA Ideal Point Distance with the USA",
                             "Log(Reserves over Current GDP)",
                             "Log(Partner Export Dependence, Current GDP)",
                             "UNGA Ideal Point Distance with the CHN"),
       reorder.coef = c(1, 15, 18, 6, 13, 2, 14, 3, 16, 4, 17, 5, c(7:12)),
       omit.coef = "factor",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol="\\wedge",
       digits = 3, 
       fontsize = "scriptsize",
       single.row = FALSE, float.pos = "!ht", sideways = TRUE,
       caption = "Results for Robustness Models Fitted to Panel Data. First five columns show logistic regression results with alternative measures of key predictors. Last two columns show conditional logistic regression results for full models. Dichotomous variable for RMB reserves in year t +  1 as the dependent variable. HAC standard errors in are parentheses for logistic regression models. Normal standard errors for conditional regression models. Year fixed-effects included in model but results omitted to facilitate presentation.", 
       label = "tb:fm_alt_measures", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)


##############################################################################
## TABLE C.3. Results from Main Logistic Regression Models Fitted to Panel Data
##############################################################################
texreg(list(io.logit.us.cs.texreg, 
            io.logit.cn.cs.texreg, 
            ti.logit.cs.texreg, 
            oc.logit.cs.texreg, 
            qpq.logit.cs.texreg, 
            full.logit.us.cs.texreg,
            full.logit.cn.cs.texreg), 
       custom.model.names = c("IO (US)", "IO (CN)", 
                              "TN", "OP", "IC", 
                              "Full (US)", "Full (CN)"),
       custom.coef.names = c("Intercept",  
                             "UNGA Ideal Point Distance with USA",
                             "BSA",
                             "Regime Type",
                             "GDP (log)",
                             "GDP per Capita (log)",
                             "Distance",
                             "Cumulative Irregular CB Turnovers",
                             "UNGA Ideal Point Distance with China",
                             "Partner Import Dependence, total imports (log)",
                             "Reserves in Months of Imports (log)",
                             "Partner Export Dependence, total exports (log)",
                             "Partner FDI Inflow Dependence"),
       reorder.coef = c(1, 2, c(9:13), c(3:8)),
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol="\\wedge",
       digits = 3, 
       fontsize = "footnotesize",
       single.row = FALSE, float.pos = "!ht", sideways = TRUE,
       caption = "Results for Logistic Regression Models Fitted to Cross-Sectional Data. Dichotomous variable for RMB reserves in year t +  1 as the dependent variable. Values for covariates are calculated as five-year averages for continuous variables. For dichotomous covariates, the covariate is assigned one for any occurrences during the sample period and zero otherwise. HAC standard errors are in parentheses. Year fixed effects included in model but results omitted to facilitate presentation.", 
       label = "tb:fm_cs", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)


##############################################################################
## FIG C.1. Comparing Logistic Regression Model Fit using ROC Curves. 
##############################################################################
# average fitted values for each model 
# io
fv.io.us <- ldply(1:length(data.mi$imputations), function(x) {
  out <- io.logit.us[[x]]$fitted.values
})
fv.io.us.mean <- apply(fv.io.us, 2, mean)
fv.io.cn <- ldply(1:length(data.mi$imputations), function(x) {
  out <- io.logit.cn[[x]]$fitted.values
})
fv.io.cn.mean <- apply(fv.io.cn, 2, mean)

# ti
fv.ti <- ldply(1:length(data.mi$imputations), function(x) {
  out <- ti.logit[[x]]$fitted.values
})
fv.ti.mean <- apply(fv.ti, 2, mean)

# oc
fv.oc <- ldply(1:length(data.mi$imputations), function(x) {
  out <- oc.logit[[x]]$fitted.values
})
fv.oc.mean <- apply(fv.oc, 2, mean)

# qpq
fv.qpq <- ldply(1:length(data.mi$imputations), function(x) {
  out <- qpq.logit[[x]]$fitted.values
})
fv.qpq.mean <- apply(fv.qpq, 2, mean)

# full ipd_us
fv.full.ipd.us <- ldply(1:length(data.mi.UN), function(x) {
  out <- full.logit.us[[x]]$fitted.values
})
fv.full.ipd.us.mean <- apply(fv.full.ipd.us, 2, mean)

# full ipd_cn
fv.full.ipd.cn <- ldply(1:length(data.mi.UN), function(x) {
  out <- full.logit.cn[[x]]$fitted.values
})
fv.full.ipd.cn.mean <- apply(fv.full.ipd.cn, 2, mean)

pdf(file = paste(FIG_DIR, "ROCplots.pdf", sep = "/"), width = 10, height = 8)
par(mfrow = c(2, 3), mai = c(0.8, 0.6, 0.6, 0.2))
# Full (US) model vs. Full (CN) model
rocplot(full.logit.us[[1]]$y, full.logit.cn[[1]]$y, fv.full.ipd.us.mean, fv.full.ipd.cn.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "Full Model (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "Full Model (CN)", pos = 4, cex  = 0.8)

# Full (US) model vs. IO (US) model
rocplot(full.logit.us[[1]]$y, io.logit.us[[1]]$y, fv.full.ipd.us.mean, fv.io.us.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "Full Model (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "IO (US) Model", pos = 4, cex  = 0.8)

# Full (US) model vs. IO (CN) model
rocplot(full.logit.us[[1]]$y, io.logit.cn[[1]]$y, fv.full.ipd.us.mean, fv.io.cn.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "Full Model (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "IO (CN) Model", pos = 4, cex  = 0.8)

# IO (US) model vs. TI model
rocplot(io.logit.us[[1]]$y, ti.logit[[1]]$y, fv.io.us.mean, fv.ti.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "IO (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "TN Model", pos = 4, cex  = 0.8)

# IO (US) model vs. OC model
rocplot(io.logit.us[[1]]$y, oc.logit[[1]]$y, fv.io.us.mean, fv.oc.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "IO (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "OP Model", pos = 4, cex  = 0.8)

# IO (US) model vs. QPQ model
rocplot(io.logit.us[[1]]$y, qpq.logit[[1]]$y, fv.io.us.mean, fv.qpq.mean)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "IO (US)", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "IC Model", pos = 4, cex  = 0.8)
dev.off()

##############################################################################
## FIG C.2. Average Ideal Point Distance with the USA and CHN, 2009-2012.
##############################################################################
# subset and summarize data
p.ipd.ave.df <- data %>%
  filter(year != 2013) %>% 
  filter(!iso3c == "USA") %>%
  select(iso3c, year, ipd_us, ipd_cn, rmbres_lead) %>%
  group_by(iso3c) %>%
  summarise(ipd_us = mean(ipd_us, na.rm = TRUE),
            ipd_cn = mean(ipd_cn, na.rm = TRUE),
            rmbres_lead = ifelse(any(rmbres_lead == 1), 1, 0))

# set parameters
shift.x <- -0.15
shift.y <- 0

# plot
p.ipd <- ggplot(p.ipd.ave.df, aes(x = ipd_us, y = ipd_cn, 
                                  group = as.factor(rmbres_lead), 
                                  color = as.factor(rmbres_lead),
                                  shape = as.factor(rmbres_lead),
                                  label = iso3c)) + 
  geom_point(alpha = 0.7, size = 6) + 
  scale_colour_manual(values = c("darkgrey", I(hsv(0/12, 7/12, 7/12)))) + 
  scale_shape_manual(values = c(1, 16)) + 
  xlab("\n UNGA Ideal Point Distance (USA)") + ylab("UNGA Ideal Point Distance (CHN)\n") + 
  scale_x_continuous(breaks = seq(0, 4.5, 0.5), limits = c(0, 4.5)) + 
  scale_y_continuous(breaks = seq(0, 4.5, 0.5), limits = c(0, 4.5)) + 
  theme_bw() + theme(axis.title.y = element_text(size = 16),
                     axis.title.x = element_text(size = 16),
                     axis.text.y = element_text(size = 14),
                     axis.text.x = element_text(size = 14),
                     legend.position = "none")

# extract country labels
label.cty <- filter(p.ipd.ave.df, rmbres_lead == 1 & !is.na(ipd_us))$iso3c

# annotate country names
for(i in 1:length(label.cty)){
  p.ipd <- p.ipd + annotate("text", x = filter(p.ipd.ave.df, iso3c == label.cty[[i]])$ipd_us + shift.x, y = filter(p.ipd.ave.df, iso3c == label.cty[[i]])$ipd_cn + shift.y, size = 3, label = label.cty[[i]])
}

pdf(file = paste(FIG_DIR, "scatter_ipd.pdf", sep = "/"), width = 10, height = 10)
p.ipd
dev.off()
